Continental Scale Hydrostratigraphy: Basin-Scale Testing of Alternative Data-Driven Approaches.

Integrated hydrological modeling is an effective method for understanding interactions between parts of the hydrologic cycle, quantifying water resources, and furthering knowledge of hydrologic processes. However, these models are dependent on robust and accurate datasets that physically represent spatial characteristics as model inputs. This study evaluates multiple data-driven approaches for estimating hydraulic conductivity and subsurface properties at the continental-scale, constructed from existing subsurface dataset components. Each subsurface configuration represents upper (unconfined) hydrogeology, lower (confined) hydrogeology, and the presence of a vertical flow barrier. Configurations are tested in two large-scale U.S. watersheds using an integrated model. Model results are compared to observed streamflow and steady state water table depth (WTD). We provide model results for a range of configurations and show that both WTD and surface water partitioning are important indicators of performance. We also show that geology data source, total subsurface depth, anisotropy, and inclusion of a vertical flow barrier are the most important considerations for subsurface configurations. While a range of configurations proved viable, we provide a recommended Selected National Configuration 1 km resolution subsurface dataset for use in distributed large-and continental-scale hydrologic modeling.

Continental Scale Hydrostratigraphy: Comparing Geologically Informed Data Products to Analytical Solutions.

This study synthesizes two different methods for estimating hydraulic conductivity (K) at large scales. We derive analytical approaches that estimate K and apply them to the contiguous United States. We then compare these analytical approaches to three-dimensional, national gridded K data products and three transmissivity (T) data products developed from publicly available sources. We evaluate these data products using multiple approaches: comparing their statistics qualitatively and quantitatively and with hydrologic model simulations. Some of these datasets were used as inputs for an integrated hydrologic model of the Upper Colorado River Basin and the comparison of the results with observations was used to further evaluate the K data products. Simulated average daily streamflow was compared to daily flow data from 10 USGS stream gages in the domain, and annually averaged simulated groundwater depths are compared to observations from nearly 2000 monitoring wells. We find streamflow predictions from analytically informed simulations to be similar in relative bias and Spearman's rho to the geologically informed simulations. R-squared values for groundwater depth predictions are close between the best performing analytically and geologically informed simulations at 0.68 and 0.70 respectively, with RMSE values under 10 m. We also show that the analytical approach derived by this study produces estimates of K that are similar in spatial distribution, standard deviation, mean value, and modeling performance to geologically-informed estimates. The results of this work are used to inform a follow-on study that tests additional data-driven approaches in multiple basins within the contiguous United States.

Water Table Depth Estimates over the Contiguous United States Using a Random Forest Model.

Water table depth (WTD) has a substantial impact on the connection between groundwater dynamics and land surface processes. Due to the scarcity of WTD observations, physically-based groundwater models are growing in their ability to map WTD at large scales; however, they are still challenged to represent simulated WTD compared to well observations. In this study, we develop a purely data-driven approach to estimating WTD at continental scale. We apply a random forest (RF) model to estimate WTD over most of the contiguous United States (CONUS) based on available WTD observations. The estimated WTD are in good agreement with well observations, with a Pearson correlation coefficient (r) of 0.96 (0.81 during testing), a Nash-Sutcliffe efficiency (NSE) of 0.93 (0.65 during testing), and a root mean square error (RMSE) of 6.87 m (15.31 m during testing). The location of each grid cell is rated as the most important feature in estimating WTD over most of the CONUS, which might be a surrogate for spatial information. In addition, the uncertainty of the RF model is quantified using quantile regression forests. High uncertainties are generally associated with locations having a shallow WTD. Our study demonstrates that the RF model can produce reasonable WTD estimates over most of the CONUS, providing an alternative to physics-based modeling for modeling large-scale freshwater resources. Since the CONUS covers many different hydrologic regimes, the RF model trained for the CONUS may be transferrable to other regions with a similar hydrologic regime and limited observations.

A hydrological simulation dataset of the Upper Colorado River Basin from 1983 to 2019.

This article presents a hydrological reconstruction of the Upper Colorado River Basin with an hourly temporal resolution, and 1-km spatial resolution from October 1982 to September 2019. The validated dataset includes a suite of hydrologic variables including streamflow, water table depth, snow water equivalent (SWE) and evapotranspiration (ET) simulated by an integrated hydrological model, ParFlow-CLM. The dataset was validated over the period with a combination of point observations and remotely sensed products. These datasets provide a long-term, natural-flow, simulation for one of the most over-allocated basins in the world.

Modeling geogenic and atmospheric nitrogen through the East River Watershed, Colorado Rocky Mountains.

There is a growing understanding of the role that bedrock weathering can play as a source of nitrogen (N) to soils, groundwater and river systems. The significance is particularly apparent in mountainous environments where weathering fluxes can be large. However, our understanding of the relative contributions of rock-derived, or geogenic, N to the total N supply of mountainous watersheds remains poorly understood. In this study, we develop the High-Altitude Nitrogen Suite of Models (HAN-SoMo), a watershed-scale ensemble of process-based models to quantify the relative sources, transformations, and sinks of geogenic and atmospheric N through a mountain watershed. Our study is based in the East River Watershed (ERW) in the Upper Colorado River Basin. The East River is a near-pristine headwater watershed underlain primarily by an N-rich Mancos Shale bedrock, enabling the timing and magnitude of geogenic and atmospheric contributions to watershed scale dissolved N-exports to be quantified. Several calibration scenarios were developed to explore equifinality using >1600 N concentration measurements from streams, groundwater, and vadose zone samples collected over the course of four years across the watershed. When accounting for recycling of N through plant litter turnover, rock weathering accounts for approximately 12% of the annual dissolved N sources to the watershed in the most probable calibration scenario (0-31% in other scenarios), and 21% (0-44% in other scenarios) when considering only "new" N sources (i.e. geogenic and atmospheric). On an annual scale, instream dissolved N elimination, plant turnover (including cattle grazing) and atmospheric deposition are the most important controls on N cycling.

Simulating Groundwater-Streamflow Connections in the Upper Colorado River Basin.

In mountain, snow driven catchments, snowmelt is supposed to be the primary contribution to river streamflows during spring. In these catchments the contribution of groundwater is not well documented because of the difficulty to monitor groundwater in such complex environment with deep aquifers. In this study we use an integrated hydrologic model to conduct numerical experiments that help quantify the effect of lateral groundwater flow on total annual and peak streamflow in predevelopment conditions. Our simulations focus on the Upper Colorado River Basin (UCRB; 2.8 × 105 km2 ) a well-documented mountain catchment for which both streamflow and water table measurements are available for several important sub-basins. For the simulated water year, our results suggest an increase in peak flow of up to 57% when lateral groundwater flow processes are included-an unexpected result for flood conditions generally assumed independent of groundwater. Additionally, inclusion of lateral groundwater flow moderately improved the model match to observations. The correlation coefficient for mean annual flows improved from 0.84 for the no lateral groundwater flow simulation to 0.98 for the lateral groundwater flow one. Spatially we see more pronounced differences between lateral and no lateral groundwater flow cases in areas of the domain with steeper topography. We also found distinct differences in the magnitude and spatial distribution of streamflow changes with and without lateral groundwater flow between Upper Colorado River Sub-basins. A sensitivity test that scaled hydraulic conductivity over two orders of magnitude was conducted for the lateral groundwater flow simulations. These results show that the impact of lateral groundwater flow is as large or larger than an order of magnitude change in hydraulic conductivity. While our results focus on the UCRB, we feel that these simulations have relevance to other headwaters systems worldwide.

Integrated Hydrologic Modeling to Untangle the Impacts of Water Management During Drought.

Over the past century, groundwater levels in California's San Joaquin Valley have dropped by more than 30 m in some areas mostly due to excessive groundwater extraction used to irrigate agricultural lands and sustain a growing population. Between 2012 and 2015, California experienced the worst drought in its recorded history, depleting surface water supplies and further exacerbating groundwater depletion in the region. Due to a lack of groundwater regulation, exact quantities of extracted groundwater in California are unknown and hard to quantify. Recent adoption of the Sustainable Groundwater Management Act has intensified efforts to identify sustainable groundwater use. However, understanding sustainable use in a highly productive agricultural system with an extremely complex surface water allocation system, variable groundwater use, and spatially extensive and diverse irrigation practices is no easy task. Using an integrated hydrologic model coupled with a land surface model, we evaluated how water management activities, specifically a suite of irrigation and groundwater pumping scenarios, impact surface water-groundwater fluxes and storage components and how those activities and the relationships between them change during drought. Results showed that groundwater pumping volume had the most significant impact on long-term water storage changes. A comparison with total water storage anomaly (TWSA) estimates from NASA's Gravity Recover and Climate Experiment (GRACE) provided some insight regarding which combinations of pumping and irrigation matched the GRACE TWSA estimates, lending credibility to these scenarios. In addition, the majority of long-term water storage changes during the recent drought occurred in groundwater storage in the deeper subsurface.

Evapotranspiration depletes groundwater under warming over the contiguous United States.

A warmer climate increases evaporative demand. However, response to warming depends on water availability. Existing earth system models represent soil moisture but simplify groundwater connections, a primary control on soil moisture. Here we apply an integrated surface-groundwater hydrologic model to evaluate the sensitivity of shallow groundwater to warming across the majority of the US. We show that as warming shifts the balance between water supply and demand, shallow groundwater storage can buffer plant water stress; but only where shallow groundwater connections are present, and not indefinitely. As warming persists, storage can be depleted and connections lost. Similarly, in the arid western US warming does not result in significant groundwater changes because this area is already largely water limited. The direct response of shallow groundwater storage to warming demonstrates the strong and early effect that low to moderate warming may have on groundwater storage and evapotranspiration.

Simulating the sensitivity of evapotranspiration and streamflow to large-scale groundwater depletion.

Groundwater pumping has caused marked aquifer storage declines over the past century. In addition to threatening the viability of groundwater-dependent economic activities, storage losses reshape the hydrologic landscape, shifting groundwater surface water exchanges and surface water availability. A more comprehensive understanding of modern groundwater-depleted systems is needed as we strive for improved simulations and more efficient water resources management. Here, we begin to address this gap by evaluating the impact of 100 years of groundwater declines across the continental United States on simulated watershed behavior. Subsurface storage losses reverberate throughout hydrologic systems, decreasing streamflow and evapotranspiration. Evapotranspiration declines are focused in water-limited periods and shallow groundwater regions. Streamflow losses are widespread and intensify along drainage networks, often occurring far from the point of groundwater abstraction. Our integrated approach illustrates the sensitivity of land surface simulations to groundwater storage levels and a path toward evaluating these connections in large-scale models.

Predicting Chronic Climate-Driven Disturbances and Their Mitigation.

Society increasingly demands the stable provision of ecosystem resources to support our population. Resource risks from climate-driven disturbances, including drought, heat, insect outbreaks, and wildfire, are growing as a chronic state of disequilibrium results from increasing temperatures and a greater frequency of extreme events. This confluence of increased demand and risk may soon reach critical thresholds. We explain here why extreme chronic disequilibrium of ecosystem function is likely to increase dramatically across the globe, creating no-analog conditions that challenge adaptation. We also present novel mechanistic theory that combines models for disturbance mortality and metabolic scaling to link size-dependent plant mortality to changes in ecosystem stocks and fluxes. Efforts must anticipate and model chronic ecosystem disequilibrium to properly prepare for resilience planning.

Connections between groundwater flow and transpiration partitioning.

Understanding freshwater fluxes at continental scales will help us better predict hydrologic response and manage our terrestrial water resources. The partitioning of evapotranspiration into bare soil evaporation and plant transpiration remains a key uncertainty in the terrestrial water balance. We used integrated hydrologic simulations that couple vegetation and land-energy processes with surface and subsurface hydrology to study transpiration partitioning at the continental scale. Both latent heat flux and partitioning are connected to water table depth, and including lateral groundwater flow in the model increases transpiration partitioning from 47 ± 13 to 62 ± 12%. This suggests that lateral groundwater flow, which is generally simplified or excluded in Earth system models, may provide a missing link for reconciling observations and global models of terrestrial water fluxes.

Using geochemical indicators to distinguish high biogeochemical activity in floodplain soils and sediments.

A better understanding of how microbial communities interact with their surroundings in physically and chemically heterogeneous subsurface environments will lead to improved quantification of biogeochemical reactions and associated nutrient cycling. This study develops a methodology to predict potential elevated rates of biogeochemical activity (microbial "hotspots") in subsurface environments by correlating microbial DNA and aspects of the community structure with the spatial distribution of geochemical indicators in subsurface sediments. Multiple linear regression models of simulated precipitation leachate, HCl and hydroxylamine extractable iron and manganese, total organic carbon (TOC), and microbial community structure were used to identify sample characteristics indicative of biogeochemical hotspots within fluvially-derived aquifer sediments and overlying soils. The method has been applied to (a) alluvial materials collected at a former uranium mill site near Rifle, Colorado and (b) relatively undisturbed floodplain deposits (soils and sediments) collected along the East River near Crested Butte, Colorado. At Rifle, 16 alluvial samples were taken from 8 sediment cores, and at the East River, 46 soil/sediment samples were collected across and perpendicular to 3 active meanders and an oxbow meander. Regression models using TOC and TOC combined with extractable iron and manganese results were determined to be the best fitting statistical models of microbial DNA (via 16S rRNA gene analysis). Fitting these models to observations in both contaminated and natural floodplain deposits, and their associated alluvial aquifers, demonstrates the broad applicability of the geochemical indicator based approach.

Metal fate and partitioning in soils under bark beetle-killed trees.

Recent mountain pine beetle infestation in the Rocky Mountains of North America has killed an unprecedented acreage of pine forest, creating an opportunity to observe an active re-equilibration in response to widespread land cover perturbation. This work investigates metal mobility in beetle-impacted forests using parallel rainwater and acid leaches to estimate solid-liquid partitioning coefficients and a complete sequential extraction procedure to determine how metals are fractionated in soils under trees experiencing different phases of mortality. Geochemical model simulations analyzed in consideration with experimental data provide additional insight into the mechanisms controlling metal complexation. Metal and base-cation mobility consistently increased in soils under beetle-attacked trees relative to soil under healthy trees. Mobility increases were more pronounced on south facing slopes and more strongly correlated to pH under attacked trees than under healthy trees. Similarly, soil moisture was significantly higher under dead trees, related to the loss of transpiration and interception. Zinc and cadmium content increased in soils under dead trees relative to living trees. Cadmium increases occurred predominantly in the exchangeable fraction, indicating increased mobilization potential. Relative increases of zinc were greatest in the organic fraction, the only fraction where increases in copper were observed. Model results reveal that increased organic complexation, not changes in pH or base cation concentrations, can explain the observed differences in metal partitioning for zinc, nickel, cadmium, and copper. Predicted concentrations would be unlikely to impair human health or plant growth at these sites; however, higher exchangeable metals under beetle-killed trees relative to healthy trees suggest a possible decline in riverine ecosystem health and water quality in areas already approaching criteria limits and drinking water standards. Impairment of water quality in important headwater streams from the increased potential for metal mobilization and storage will continue to change as beetle-killed trees decompose and forests begin to recover.

The effects of physical and geochemical heterogeneities on hydro-geochemical transport and effective reaction rates.

The role of coupled physical and geochemical heterogeneities in hydro-geochemical transport is investigated by simulating three-dimensional transport in a heterogeneous system with kinetic mineral reactions. Ensembles of 100 physically heterogeneous realizations were simulated for three geochemical conditions: 1) spatially homogeneous reactive mineral surface area, 2) reactive surface area positively correlated to hydraulic heterogeneity, and 3) reactive surface area negatively correlated to hydraulic heterogeneity. Groundwater chemistry and the corresponding effective reaction rates were calculated at three transverse planes to quantify differences in plume evolution due to heterogeneity in mineral reaction rates and solute residence time (τ). The model is based on a hypothetical CO2 intrusion into groundwater from a carbon capture utilization and storage (CCUS) operation where CO2 dissolution and formation of carbonic acid created geochemical dis-equilibrium between fluids and the mineral galena that resulted in increased aqueous lead (Pb(2+)) concentrations. Calcite dissolution buffered the pH change and created conditions of galena oversaturation, which then reduced lead concentrations along the flow path. Near the leak kinetic geochemical reactions control the release of solutes into the fluid, but further along the flow path mineral solubility controls solute concentrations. Simulation results demonstrate the impact of heterogeneous distribution of geochemical reactive surface area in coordination with physical heterogeneity on the effective reaction rate (Krxn,eff) and Pb(2+) concentrations within the plume. Dissimilarities between ensemble Pb(2+) concentration and Krxn,eff are attributed to how geochemical heterogeneity affects the time (τeq) and therefore advection distance (Leq) required for the system to re-establish geochemical equilibrium. Only after geochemical equilibrium is re-established, Krxn,eff and Pb(2+) concentrations are the same for all three geochemical conditions. Correlation between reactive surface area and hydraulic conductivity, either positive or negative, results in variation in τeq and Leq.

Human health risk assessment of CO2 leakage into overlying aquifers using a stochastic, geochemical reactive transport approach.

Increased human health risk associated with groundwater contamination from potential carbon dioxide (CO2) leakage into a potable aquifer is predicted by conducting a joint uncertainty and variability (JUV) risk assessment. The approach presented here explicitly incorporates heterogeneous flow and geochemical reactive transport in an efficient manner and is used to evaluate how differences in representation of subsurface physical heterogeneity and geochemical reactions change the calculated risk for the same hypothetical aquifer scenario where a CO2 leak induces increased lead (Pb(2+)) concentrations through dissolution of galena (PbS). A nested Monte Carlo approach was used to take Pb(2+) concentrations at a well from an ensemble of numerical reactive transport simulations (uncertainty) and sample within a population of potentially exposed individuals (variability) to calculate risk as a function of both uncertainty and variability. Pb(2+) concentrations at the well were determined with numerical reactive transport simulation ensembles using a streamline technique in a heterogeneous 3D aquifer. Three ensembles with variances of log hydraulic conductivity (σ(2)lnK) of 1, 3.61, and 16 were simulated. Under the conditions simulated, calculated risk is shown to be a function of the strength of subsurface heterogeneity, σ(2)lnK and the choice between calculating Pb(2+) concentrations in groundwater using equilibrium with galena and kinetic mineral reaction rates. Calculated risk increased with an increase in σ(2)lnK of 1 to 3.61, but decreased when σ(2)lnK was increased from 3.61 to 16 for all but the highest percentiles of uncertainty. Using a Pb(2+) concentration in equilibrium with galena under CO2 leakage conditions (PCO2 = 30 bar) resulted in lower estimated risk than the simulations where Pb(2+) concentrations were calculated using kinetic mass transfer reaction rates for galena dissolution and precipitation. This study highlights the importance of understanding both hydrologic and geochemical conditions when numerical simulations are used to perform quantitative risk calculations.

A new perspective on human health risk assessment: development of a time dependent methodology and the effect of varying exposure durations.

We present a new Time Dependent Risk Assessment (TDRA) that stochastically considers how joint uncertainty and inter-individual variability (JUV) associated with human health risk change as a function of time. In contrast to traditional, time independent assessments of risk, this new formulation relays information on when the risk occurs, how long the duration of risk is, and how risk changes with time. Because the true exposure duration (ED) is often uncertain in a risk assessment, we also investigate how varying the magnitude of fixed size durations (ranging between 5 and 70 years) of this parameter affects the distribution of risk in both the time independent and dependent methodologies. To illustrate this new formulation and to investigate these mechanisms for sensitivity, an example of arsenic contaminated groundwater is used in conjunction with two scenarios of different environmental concentration signals resulting from rate dependencies in geochemical reactions. Cancer risk is computed and compared using environmental concentration ensembles modeled with sorption as 1) a linear equilibrium assumption (LEA) and 2) first order kinetics (Kin). Results show that the information attained in the new time dependent methodology reveals how the uncertainty in other time-dependent processes in the risk assessment may influence the uncertainty in risk. We also show that individual susceptibility also affects how risk changes in time, information that would otherwise be lost in the traditional, time independent methodology. These results are especially pertinent for forecasting risk in time, and for risk managers who are assessing the uncertainty of risk.

Kinetic metal release from competing processes in aquifers.

Understanding groundwater time scales wherein kinetic metal-desorption and mineral-dissolution are important mechanisms is essential for realistic modeling of metal release. In this study, release rate constants were compiled and the Damköhler number was applied to calculate residence times where kinetic formulations are relevant. Desorption rate constants were compiled for arsenic, barium, cadmium, copper, lead, mercury, nickel, and zinc, and span 6 orders of magnitude, while mineral-dissolution rate constants compiled for calcite, kaolinite, smectite, anorthite, albite, K-feldspar, muscovite, quartz, goethite, and galena ranged over 13 orders of magnitude. This Damköhler analysis demonstrated that metal-desorption kinetics are potentially influential at residence times up to about two years, depending on the metal and groundwater conditions. Kinetic mineral-dissolution should be considered for nearly all residence times relevant to groundwater modeling, provided the rate, solubility, and availability of the mineral generates a non-negligible concentration. Geochemical models of competitive desorption and dissolution for an illustrative metal demonstrate total metal concentrations may be sensitive to dissolution rate variations despite the predominance of release from desorption. Ultimately, this analysis provides constraints on relevant processes for incorporation into transport models.